Corneal dynamic model algorithm and system using the same

ABSTRACT

A corneal dynamic model algorithm and system using the same is provided, comprising an ophthalmotonometer and a calculating unit. The calculating unit comprises the algorithm for calculating the corneal dynamic model. The algorithm comprises: reading an actual corneal deformation from the ophthalmotonometer; substituting a young&#39;s modulus initial value and a damping coefficient initial value into a mathematical equation for getting a calculated corneal deformation; determining whether error amount between the calculated deformation and the actual deformation is a minimum or not; and getting a young&#39;s modulus and a damping coefficient.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a non-provisional application which claims priority of U.S. Provisional Application No. 62/134,438 filed on Mar. 17, 2015 under 35 U.S.C. § 119(e), the entire contents of all of which are hereby incorporated by reference.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to a corneal dynamic model algorithm and system using the same. Specially, it is an instantaneous and non-invasive corneal dynamic model algorithm and system using the same.

2. Description of the Prior Art

As technology advances, it also drives the speedy development of computers, televisions, mobile phones and other digital merchandises. The long hours of using eyes on surfing internet, watching TV, operating the mobile phone cause eyes dry, photophobia, tearing, fatigue, and even with conditions that are accompanied with symptoms of headache, dizziness, nausea, shoulder and neck pain, blurred vision. As a result, eye fatigue has become a common civilization disease of people in the modern society.

Overuse of eyes can result in abnormal function of focus adjustment, loss of efficacy in adjusting the focus length precisely, blurred vision, aggravating myopia, xerophthalmia, or eye diseases such as induced glaucoma and retinal lesion.

Therefore, in addition to a timely rest for the eyes, regular examination of the eyes is particularly important. The advanced intra-ocular pressure (IOP) detector is the Corvis® ST of Oculus in Germany. The Corvis® ST is a non-invasive ophthalmotonometer. It uses the high-speed photographic technology to snap 140 corneal sectional images within 31 ms, and calculate and analyze data, including the degree of applanation and rebound when a specific amount of air pressure is blown onto the eye and the speed, in order to evaluate biomechanical characteristics of the cornea.

So far the Corvis® ST in medical science only applies to detect IOP, cornealectasia, and collagen crosslinking. However, it cannot detect instantaneously the Young's modulus of cornea.

SUMMARY OF THE INVENTION

In view of the above problems, in one aspect, the present invention provides a dynamic model for analyzing the mechanical properties instantaneously when receiving the images from the ophthalmotonometer.

In one embodiment, the present invention provides a corneal dynamic model algorithm, comprising: (S1) reading an actual corneal deformation from an ophthalmotonometer; (S2) substituting a young's modulus initial value and a damping coefficient initial value into a mathematical equation for getting a calculated corneal deformation; (S3) determining whether error amount between the calculated deformation and the actual deformation is a minimum or not; and (S4) getting a young's modulus and a damping coefficient.

Wherein if the determining in the step (S3) is not, then back to the step (S2) and iterate in optimization until the error amount is a minimum.

In this embodiment, the mathematical equation is

${w\left( {\theta,{\phi;t}} \right)} = {\sum\limits_{l = 1}^{\infty}\;{{w_{l\; 0}\left( {\theta,{\phi;t}} \right)}.}}$

In this embodiment, the minimum is 10⁻³.

In another aspect, the present invention provides a corneal dynamic model system, comprising an ophthalmotonometer and a calculating unit. The calculating unit comprises the algorithm for calculating the corneal dynamic model. The algorithm comprises: (S1) reading an actual corneal deformation from the ophthalmotonometer; (S2) substituting a young's modulus initial value and a damping coefficient initial value into a mathematical equation for getting a calculated corneal deformation; (S3) determining whether error amount between the calculated deformation and the actual deformation is a minimum or not; and (S4) getting a young's modulus and a damping coefficient.

Wherein if the determining in the step (S3) is not, then back to the step (S2) and iterate in optimization until the error amount is a minimum.

In this embodiment, the mathematical equation is

${w\left( {\theta,{\phi;t}} \right)} = {\sum\limits_{l = 1}^{\infty}\;{{w_{l\; 0}\left( {\theta,{\phi;t}} \right)}.}}$

In this embodiment, the minimum is 10⁻³.

The following description and the drawings set forth certain illustrative advantages and spirits of the specification.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A is a schematic diagram of an eyeball.

FIG. 1B is a schematic diagram of an eyeball dynamic model in one embodiment.

FIG. 2A-FIG. 2F are schematic diagrams of the first six model shapes of the eyeball model in one embodiment.

FIG. 3 is a schematic diagram of a corneal dynamic model system in one embodiment.

FIG. 4A is a schematic diagram of air pressure Gaussian distribution.

FIG. 4B is a schematic diagram of time series of the air pressure.

FIG. 5A is a schematic diagram of time series displacement at the top of cornea.

FIG. 5B is a schematic diagram of the first six modal time series displacement at the top of cornea.

FIG. 6 is a flowchart of a corneal dynamic model algorithm in one embodiment.

FIG. 7 is corneal images taken from one patient by the ophthalmotonometer.

FIG. 8A is a schematic diagram of one patient's corneal displacement time series from the ophthalmotonometer and our model.

FIG. 8B is a schematic diagram of the damping ratio affections.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

The following will disclose a plurality of embodiment with texts and drawings. For a clear illustration, many details in practice will be described in the following description. However, it should be understand that the details of the practical application of these are not to limit the present invention. In addition, for the sake of simplifying the drawings, some of the conventional structures and components in the drawings will be depicted in a simplified schematic manner.

In order to implement this invention, we simply the eyeball structure. There are eight assumptions setup through this invention. As shown in FIG. 1A, the eye is made up of the three layers and the three transparent structures. The outermost layer is composed of the cornea and sclera, the middle layer is consists of the choroid, ciliary body (not indicated), and iris, and the innermost is the retina. The eight basic assumptions are made: (i) the eyeball is a perfectly spherical diaphragm, and fluid is completely filled inside the eyeball; (ii) the thickness and material of the diaphragm are uniform everywhere. Here we take the corneal properties as the material of diaphragm; (iii) the amplitude of oscillation is small compared with the characteristic length of the eyeball; (iv) the fluid inside the eyeball is incompressible water; (v) the iris is neglected because of its relatively small volume; (vi) the space of lens is regarded as fluid because of its density is closed to fluid; (vii) the air impulsion applies on the cornea with a Gaussian distribution; (viii) boundary conditions at the surface outside the eyeball is free movement.

As shown in FIG. 1B, this model simplifies the eyeball into the three parts: the air outside the eyeball (the oblique and rectus muscles are not concerned), the diaphragm of the eyeball (including cornea, sclera and retina), and the fluid inside the eyeball (including vitreous humor, lens, and aqueous humor.) Here, we have three steps to derive the modal frequencies of the eyeball.

Firstly, the spherical surface harmonics and radiation functions help develop the pressure functions outside and inside the eyeball. The Euler's equation for the ambient air and the fluid inside the eyeball is

$\begin{matrix} \left\{ {\begin{matrix} {{\frac{\partial v}{\partial t} + {v \cdot {\nabla v}}} = {{- \frac{1}{\rho}}{\nabla p}}} \\ {{\nabla{\cdot v}} = 0} \end{matrix}.} \right. & \left( {{eq}.\mspace{14mu} 1} \right) \end{matrix}$ Wherein v_(i) is the covariant component of velocity fields, ρ is the density, p is the pressure, and we neglect the fluid viscosity. In the spherical coordinates system (r, θ, φ, the linearized hydrodynamic equation with velocity in r direction is

$\begin{matrix} {\frac{\partial v_{r}}{\partial t} = {{- \frac{1}{\rho}}{\frac{\partial p}{\partial r}.}}} & \left( {{eq}.\mspace{14mu} 2} \right) \end{matrix}$ The solution of p⁺ inside the eyeball is given by

$\begin{matrix} {p^{+} = {\sum\limits_{l = 0}^{\infty}\;{\sum\limits_{m = l}^{- l}\;{i\;\rho^{+}R^{2}\omega^{2}A_{lm}^{+}{Y_{lm}\left( {\theta,\phi} \right)}\left( \frac{r}{R} \right)^{l}{e^{i\;\omega\; t}.}}}}} & \left( {{eq}.\mspace{14mu} 3} \right) \end{matrix}$ With undetermined constants A_(lm) ⁺ and where the angular frequency ω, R is the radius of eyeball, and Y_(lm) are the spherical surface harmonics.

Note that the superscript ‘+’ denotes the domain inside the eyeball. The Fourier transform of functions with respective to time t is given by {circumflex over (x)}(ω)=∫_(−∞) ^(∞)x(t)e^(iωt)dt. So that the solution of p⁻ outside the eyeball is

$\begin{matrix} {p^{-} = {\sum\limits_{l = 0}^{\infty}\;{\sum\limits_{m = l}^{- l}\;{i\;\rho^{-}R^{2}\omega^{2}A_{lm}^{-}{Y_{lm}\left( {\theta,\phi} \right)}\left( \frac{r}{R} \right)^{{- l} - 1}{e^{i\;\omega\; t}.}}}}} & \left( {{eq}.\mspace{14mu} 4} \right) \end{matrix}$ Here the air and fluid are assumed to be incompressible, so the zeroth harmonic Y₀₀(θ, ϕ) is excluded. The velocity component on the surface, i.e. r=R, becomes

$\begin{matrix} \begin{matrix} {{V_{R}\left( {r,\theta,{\phi;t}} \right)} = {\sum\limits_{l = 0}^{\infty}\;{\sum\limits_{m = l}^{- l}{{- \rho^{-}}R^{2}\omega\; A_{lm}^{+}{Y_{lm}\left( {\theta,\phi} \right)}e^{i\;\omega\; t}}}}} \\ {= {\sum\limits_{l = 0}^{\infty}\;{\sum\limits_{m = l}^{- l}{{- \rho^{-}}R^{2}\omega\; A_{lm}^{-}{Y_{lm}\left( {\theta,\phi} \right)}{e^{i\;\omega\; t}.}}}}} \end{matrix} & \left( {{eq}.\mspace{14mu} 5} \right) \end{matrix}$ Note that the series expansion of both air and fluid are substantially setup on the orthogonal and complete sets.

Secondly, the modified Helmholtz function helps derive the displacement function on the stretched diaphragm. The stretched diaphragm with tension and stiffness has a potential energy of bending which is determined by the tension of the diaphragm, fluid pressure, and the elastic constants of the material. The intraocular pressure let the diaphragm be stretched tightly, and the quantity of tension is 10-20% of the stiffness of the material of the diaphragm. Because the added tension increases the effective modulus, the added term T to the Young's modulus occurs. Hence the modified Helmholtz equation satisfied for points at the surface is

$\begin{matrix} {{{\frac{\left( {E + T} \right) \cdot t^{3}}{12\left( {1 - v^{2}} \right)}{\nabla^{4}w}} - {{T \cdot t}{\nabla^{2}w}} - \lbrack p\rbrack + {{\rho \cdot t}\frac{\partial^{2}w}{\partial t^{2}}}} = 0.} & \left( {{eq}.\mspace{14mu} 6} \right) \end{matrix}$ Wherein E is the Young's modulus, T is the tension on the diaphragm, t is the thickness of the diaphragm, v is the Poisson's ratio, ρ is the density of the diaphragm, and [p] is the pressure difference between the two sides of the diaphragm, i.e. [p]=p⁺−p⁻. For the case of a spherical diaphragm, the Laplacian of w on sphere with (r, θ, φ) are

$\begin{matrix} {\nabla^{2}{= {{\frac{1}{r^{2}}\left\lbrack {{\frac{1}{\sin\;\theta}\frac{\partial\;}{\partial\theta}\left( {\sin\;\theta\frac{\partial\;}{\partial\theta}} \right)} + {\frac{1}{\sin^{2}\theta}\frac{\partial^{2}}{\partial\phi^{2}}}} \right\rbrack}.}}} & \left( {{eq}.\mspace{14mu} 7} \right) \end{matrix}$ We assume the simple harmonic motion is the solution of this equation,

$\begin{matrix} \begin{matrix} {{w\left( {\theta,{\phi;t}} \right)} = {\sum\limits_{l = 0}^{\infty}\;{\sum\limits_{m = l}^{- l}{w_{lm}\left( {\theta,{\phi;t}} \right)}}}} \\ {= {\sum\limits_{l = 0}^{\infty}\;{\sum\limits_{m = l}^{- l}{{- i}\; B_{lm}{{Y_{lm}\left( {\theta,\phi} \right)} \cdot {e^{i\;\omega\; t}.}}}}}} \end{matrix} & \left( {{eq}.\mspace{14mu} 8} \right) \end{matrix}$ Wherein B_(lm) are the undetermined constants. Note that the series expansion of the diaphragm is substantially setup on the orthogonal and complete bases.

Thirdly, the continuities on the interfaces help to develop the characteristic function. There are two radial conditions of velocities ensure the continuous on the interfaces between the air-diaphragm and the diaphragm-liquid: v _(R) ⁺ ={dot over (w)}  (eq. 9) v _(R) ⁻ ={dot over (w)}  (eq. 10). Substituting Eqs. (4), (5), and (8) into Eqs. (9) and (10) leads the relationships with three unknowns A_(lm) ⁺, A_(lm) ⁻, and B_(lm): lRA _(lm) ⁺ +B _(lm)=0  (eq. 11) (l+1)RA _(lm) ⁻ −B _(lm)=0  (eq. 12). Moreover, using Eqs. (3) and (4), replacing the coefficients of A_(lm) ⁺ and A_(lm) ⁻ by B_(lm), and modifying the pressure difference term [p] at the diaphragm lead

$\begin{matrix} {\lbrack p\rbrack_{lm} = {{- i}\; R\;\omega^{2}B_{lm}{Y_{lm}\left( {\theta,\phi} \right)}{{e^{i\;\omega\; t}\left( {\frac{\rho^{+}}{l} + \frac{\rho^{-}}{l + 1}} \right)}.}}} & \left( {{eq}.\mspace{14mu} 13} \right) \end{matrix}$ Furthermore, the second differentiation of the displacement from Eq. (8) is

$\begin{matrix} {\frac{\partial^{2}w_{lm}}{\partial t^{2}} = {i\;\omega^{2}B_{lm}{Y_{lm}\left( {\theta,\phi} \right)}{e^{i\;\omega\; t}.}}} & \left( {{eq}.\mspace{14mu} 14} \right) \end{matrix}$ So that the pressure difference can be rewritten to

$\begin{matrix} {\lbrack p\rbrack_{lm} = {{- {R\left( {\frac{\rho^{+}}{l} + \frac{\rho^{-}}{l + 1}} \right)}}{\frac{\partial^{2}w_{lm}}{\partial t^{2}}.}}} & \left( {{eq}.\mspace{14mu} 15} \right) \end{matrix}$

Substituting Eq. (15) into Eq. (6) leads the new modal motion equation

$\begin{matrix} {{\left\lbrack {{\frac{\left( {E + T} \right) \cdot t^{3}}{12\left( {1 - v^{2}} \right)}\frac{{l^{2}\left( {l + 1} \right)}^{2}}{R^{4}}} + \frac{{l\left( {l + 1} \right)}{T \cdot t}}{R^{2}}} \right\rbrack w_{lm}} + {\quad{{\left\lbrack {{\rho \cdot t} + {R\left( {\frac{\rho^{+}}{l} + \frac{\rho^{-}}{l + 1}} \right)}} \right\rbrack\frac{\partial^{2}w_{lm}}{\partial t^{2}}} = 0.}}} & \left( {{eq}.\mspace{14mu} 16} \right) \end{matrix}$ Wherein the relation of R²∇²w_(lm)=−l(l+1)w_(lm) is used. Then it is obviously the natural frequency term becomes

$\begin{matrix} {\omega_{l}^{2} = {\frac{{\frac{\left( {E + T} \right) \cdot t^{3}}{12\left( {1 - v^{2}} \right)}\frac{{l^{3}\left( {l + 1} \right)}^{3}}{R^{4}}} + \frac{T \cdot t \cdot {l^{2}\left( {l + 1} \right)}^{2}}{R^{2}}}{{\left\lbrack {{\left( {l + 1} \right)\rho^{+}} + {l \cdot p^{-}}} \right\rbrack R} + {{l\left( {l + 1} \right)}{\rho \cdot t}}}.}} & \left( {{eq}.\mspace{14mu} 17} \right) \end{matrix}$ Thus the modal frequency is f_(l)=ω_(l)/2π. Moreover, Eq. (17) can be rewritten into

$\begin{matrix} {{{M_{l}\frac{\partial^{2}w_{lm}}{\partial t^{2}}} + {K_{l}w_{lm}}} = 0.} & \left( {{eq}.\mspace{14mu} 18} \right) \end{matrix}$ Wherein the modal masses are

${M_{l} = {{\rho \cdot t} + {R\left( {\frac{\rho^{+}}{l} + \frac{\rho^{-}}{l + 1}} \right)}}},$ and modal spring constants are

$\begin{matrix} {K_{l} = {{\frac{\left( {E + T} \right) \cdot t^{3}}{12\left( {1 - v^{2}} \right)}\frac{{l^{2}\left( {l + 1} \right)}^{2}}{R^{4}}} + {\frac{{l\left( {l + 1} \right)}{T \cdot t}}{R^{2}}.}}} & \; \end{matrix}$

In the previous discussion, the modal equation of motion is given. The forced vibration also could be extended in modal analysis. While the external force term F_(lm) the damping term C_(l){dot over (w)}_(lm) are added in Eq. (18), the equation of motion of the system becomes

$\begin{matrix} {{{M_{l}\frac{\partial^{2}w_{lm}}{\partial t^{2}}} + {C_{l}\frac{\partial w_{lm}}{\partial t}} + {K_{l}w_{lm}}} = {F_{lm}.}} & \left( {{eq}.\mspace{14mu} 19} \right) \end{matrix}$ Note that the damping comes from the viscosity of water and the structural diaphragm. But the certain values of the dampers are difficult to be obtained. So we just put damper coefficients as the undetermined coefficients, and then we will use numerical technique to estimate them.

Here the series expansion of the external force is defined by

$\begin{matrix} {{F\left( {\theta,{\phi;t}} \right)} = {\sum\limits_{l = 0}^{\infty}\;{\sum\limits_{m = {- l}}^{l}\;{F_{lm}{Y_{lm}\left( {\theta,\phi} \right)}{e^{i\;\omega\; t}.}}}}} & \left( {{eq}.\mspace{14mu} 20} \right) \end{matrix}$ Wherein

$\begin{matrix} {F_{lm} = {\frac{\int_{0}^{2\pi}{\int_{0}^{\pi}{{F\left( {\theta,{\phi;t}} \right)}{Y_{lm}\left( {\theta,\phi} \right)}\sin\;\theta\ d\;\theta\ d\;\phi}}}{\int_{0}^{2\pi}{\int_{0}^{\pi}{{Y_{lm}^{2}\left( {\theta,\phi} \right)}\sin\;\theta\ d\;\theta\ d\;\phi}}}.}} & \left( {{eq}.\mspace{14mu} 21} \right) \end{matrix}$ Then solution of Eq. (19) is

$\begin{matrix} {{{w_{lm}\left( {\theta,{\phi;t}} \right)} = {\frac{1}{M_{l}}e^{{- \xi_{l}}\omega_{Dl}t}{\int_{0}^{t}{\frac{\sin\;{\omega_{Dl}\left( {t - \tau} \right)}}{\omega_{Dl}}F_{lm}\ d\;{\tau \cdot {Y_{lm}\left( {\theta,\phi} \right)}}}}}},} & \left( {{eq}.\mspace{14mu} 22} \right) \end{matrix}$ in which the damping ratio is defined by ξ_(l)=C_(l)/M_(l) and ω_(Dl)=ω√{square root over (1−ξ²)}.

Moreover, the external force is a symmetric force about the θ=0 axis, namely F(θ, ϕ′, t)=g(θ)h(t). That means m=0, and the surface harmonics could be simplified by Y_(l0)=P_(l)(cos θ). Then F _(l0)=(l+½)∫₀ ^(π) g(θ)P _(l)(θ)sin θdθ·h(t)  (eq. 23).

The integration, ∫₀ ^(π)[P_(l)(cos θ)]² sin dθ=2/(2l+1), helps to solve the modal displacements

$\begin{matrix} {{{w_{l\; 0}\left( {\theta,{\phi;t}} \right)} = {\frac{\left( {{2l} + 1} \right)}{2M_{l}\omega_{Dl}}{{P_{l}\left( {\cos\;\theta} \right)} \cdot {\int_{0}^{\pi}{{g(\alpha)}{P_{l}(\alpha)}\sin\;\alpha\ d\;{\alpha \cdot {\int_{0}^{t}{e^{{- {\xi\omega}_{Dl}}t}\sin\;{\omega_{Dl}\left( {t - \tau} \right)}{h(t)}\ d\;\tau}}}}}}}},} & \left( {{eq}.\mspace{14mu} 24} \right) \end{matrix}$ and the total displacements at the diaphragm become

$\begin{matrix} {{w\left( {\theta,{\phi;t}} \right)} = {\sum\limits_{l = 1}^{\infty}\;{{w_{l\; 0}\left( {\theta,{\phi;t}} \right)}.}}} & \left( {{eq}.\mspace{14mu} 25} \right) \end{matrix}$

It is noted that on the solution of the displacement at the diaphragm, there are two terms which should be solved before the numerical calculation: the modal shape functions of the diaphragm, the modal external force in spatial distribution and in time series. For the modal shape functions of the diaphragm, the surface harmonics Y_(lm)(θ,ϕ) are the main terms of Eq. (8). Since the external force is axial symmetry and m=0, the modal shapes [Y_(l0)=P_(l)(cos θ)] could be obtained from the definition of the Legendre polynomial function.

One embodiment in this application, please refer to FIG. 2A˜FIG. 2F and FIG. 3. A corneal dynamic model system 11 preferably comprises an ophthalmotonometer 12 and a calculating unit 13. In this embodiment, the ophthalmotonometer 12 is Corvis® ST produced by Oculus, but not limited thereto. This ophthalmotonometer uses the high-speed photographic technology to snap 140 corneal sectional images within 31 ms, and calculate and analyze data, including the degree of applanation and rebound when a specific amount of air pressure is blown onto the eye and the speed, in order to evaluate biomechanical characteristics of the cornea. However, other ophthalmotonometers can be used in other embodiments. The calculating unit 13 can be a processor or other members with computing capability.

By directing an external force (e.g. air-puff) to the eye 21 with the ophthalmotonometer 12, the ophthalmotonometer 12 could obtain an actual corneal deformation Dt. Wherein when the eyeball is directed by the ophthalmotonometer 12, the first six model shapes of the eyeball model are shown in FIG. 2A˜FIG. 2F. The red dots mark the nodal points and are regarded as the hinged conditions at the surface of the diaphragm.

In addition, as shown in FIG. 4A, we assume the pressure distribution deduced from the ophthalmotonometer 12 is a three-dimensional Gaussian distribution. The width of the Gaussian distribution is by a parameter μ. Besides, the time series of the air pressure measured by the experiment is shown in FIG. 4B.

Table I gives the material properties of the eyeball.

TABLE I Properties Symbols Values Radius R 12 mm Thickness t 540 μm Poisson's ratio ν 0.49 Density of air ρ⁻ 1.204 kg/m³ Density of fluid inside eyeball ρ⁺ 1000 kg/m³ Density of sclera ρ 1444 kg/m³ Moreover, according to Eq. 18, the top 10 modal frequencies are listed on Table II.

TABLE II Mode No. Frequency [Hz] 1 14.363 2 33.0 3 55.657 4 82.109 5 112.336 6 146.437 7 184.576 8 226.951 9 273.771 10 325.246 Wherein the Young's modulus is assumed 0.48034 MPa and the IOP is 14.5 mmHg.

In this case, Eq. 25 helps to obtain the displacements at the point as shown in FIG. 5A. The time series of the first six modes, denoting the eigen functions of Eq. 24 by the following equation 26, as shown in FIG. 5B. The eq. 26 is w _(l)(0,0;T)=P _(l)(1)·∫₀ ^(π) g(α)P _(l)(α)sin αdα·∫ ₀ ^(t) e ^(−ξω) ^(Dl) ^(t) sin ω_(l)(t−τ)h(t)dτ  (eq. 26).

It is noted that, before substituting the above parameters into the mathematical equation (i.e. Eq. 25), we need to assume initial values respectively for the Young's modulus and the damping coefficients. Wherein the assumed initial value can be a value based on the past experience or a value provided by the research work.

In above discussion, it is obviously that both of the Young's modulus and the damping coefficients affect the distribution shape of the corneal deformation. Here we compared the simulation displacements w(0,0;t) by applying Eq. (25) with these measured from the corneal displacement. Then we applied the optimization program (e.g. Matlab) to estimate the modal spring constants, K_(l), and the modal damping coefficients, C_(l).

Finding the error amount during the optimization program. That is, determine whether error amount between the calculated deformation and the actual deformation is a minimum or not. When the error amount is minimum, the Young's modulus and the damping coefficients can be obtained. It is noted that, when the error amount is determined not to be the minimal value, re-assume initial values respectively for the Young's modulus and the damping coefficients to calculate iteratively until the error value becomes the minimum value. In this embodiment, the minimum preferably is 10⁻³, but not limited thereto.

As set forth above, the whole algorithm is shown in FIG. 6. The algorithm preferably is applied to the corneal dynamic model system above mentioned. The algorithm comprising: (S1) reading an actual corneal deformation from an ophthalmotonometer; (S2) substituting a young's modulus initial value and a damping coefficient initial value into a mathematical equation for getting a calculated corneal deformation; (S3) determining whether error amount between the calculated deformation and the actual deformation is a minimum or not; and (S4) getting a young's modulus and a damping coefficient. It is noted that, if step (S3) is determined not to be the minimal value, then return to step (S2) to calculate iteratively until the error value becomes the minimum value.

Note that this algorithm can be written/stored in the calculating unit 13 of the corneal dynamic model system by software or firmware. The calculating unit 13 could connect externally with the ophthalmotonometer. However, it could also be set in the ophthalmotonometer.

In practical application, there were 25 patients joined this study, and their information is listed on Table III.

TABLE III Estimated Patient Thickness IOP Young's Damping No Age [μm] [mmHg] modulus [MPa] coefficients 1 74 742 11 0.1178 0.0498 2 23 545 16 0.1102 0.0471 3 40 565 13 0.0103 0.0372 4 26 459 12.5 0.6644 0.0437 5 28 558 14 0.1122 0.0446 6 36 585 14.5 0.1204 0.0594 7 50 578 13 0.1411 0.1 8 26 596 14.5 0.1074 0.0355 9 39 564 15 0.1114 0.0144 10 38 525 12 0.205 0.0346 11 33 501 13 0.1301 0.083 12 27 565 14 0.7325 0.133 13 26 529 15 0.2909 0.0336 14 49 518 16 0.1347 0.035 15 24 577 16.5 0.5288 0.033 16 35 547 15 0.2651 0.0343 17 27 651 17 0.1022 0.0356 18 51 584 15.5 0.026 0.033 19 25 506 13 0.5701 0.0329 20 42 454 11.5 0.6149 0.0329 21 22 544 14 1.2427 0.0323 22 36 497 16 1.0056 0.0549 23 22 571 12.5 0.2686 0.0305 24 39 476 12 0.0739 0.0342 25 26 555 13 0.2708 0.0348 Average 34.56 551.68 13.98 0.3183 0.0456 In this study, parts of the material properties were assumed and shown on Table I, except the thicknesses. The values IOPs and the thicknesses could be measured from the ophthalmotonometer 12. As shown in FIG. 7, the displacements of the cornea of one patient at time 0.000, 9.471, 17.325, and 21.252 msec. Besides, the vibration information of cornea can be obtained by technique of image modeling and analysis from images before and during the air-puff.

FIG. 8A shows one patient's corneal displacement time series measured from the ophthalmotonometer 12, and the same time series obtained from our model are also shown. The two curves are not perfectly close, because there are too many assumptions in the model. Besides, how to setup the boundary conditions is a big issue. From the calculation experience, the Young's modulus affects almost the first half oscillation of the duration; the damping ratio affects the dispersion in the last half of the duration as shown in FIG. 8B.

The above mentioned Table III shows the results of Young's moduli and the damping coefficients among the 25 patients. The average Young's modulus is about 0.3183 Mpa and the average damping coefficient is about 0.0456. These values satisfy the acceptable ranges 0.1˜10 MPa from the tests of the human eyes (in vitro) as listed on Table IV.

TABLE IV Young's moduli No. Researchers [Mpa] Remark Years References 1 Woo, SL-Y., et al 1.8-8.1 Human cadaver 1972 FIG. 6 in [19] eyes, with I-3 days post mortem 2 Jue, Bill, and 5 Human cadaver eyes, 1986 [20] David M. aged between 45 and Maurice 90 year, obtained within 24 hour of death. 3 Wang, Hsichun, 4.2-30  Six human cadaver 1996 Table 1 et al. eye, aged 75-82 years, in [21] immersed in the 15% solution of saline, the other six ones in the dextran, aged 78-84 years. 4 Uchio, Eiichi, et  50-100 Human cadaver eye 1999 FIG. 1(A) al. enucleated between in [22] 8-24 hour and stored at 4° C., experiment conducting withing 48 hour 5 Katherine D.  6.89-110.32 Human cadaver eye 2003 Table 8 Voorhies Packaged in the saline in [23] gauze 6 Elsheikh, 0.16-0.81 36 human cadaver 2007 [24] Ahmed, Defu eye, immersed in the Wang, and storage medium David Pye Optisol GS at 4° C. within 12 hours of death and tested within 14 days. 7 Hamilton, 0.13-0.43 100 living eyes, aged 2008 [25] Kirsten E., and 18 between 30 years. David C. Pye The Orssengo-Pye algorithm is used. 8 Elsheikh, 0.15-1.15 Forty nine Human 2008 FIG. 6 in [26] Ahmed, Daad cadaver eyes Alhasso, and preserved in the Eusol Paolo Rama C, for a maximum of 12 hour of death, used in 5-7 day of preservation 9 Cartwright, 0.27-0.52 Human cadaver 2011 [27] Nathaniel E. eyes, aged between 24 Knox, John R. and 102 years, with Tyrer, and John 5-27 days post Marshall mortem, immersed in Eagle's minimal essential medium

Compared with the prior art, this application provides real-time dynamic model to estimate instantaneously the eyeball properties during the IOP measuring by the ophthalmotonometer.

Although the preferred embodiments of the present invention have been described herein, the above description is merely illustrative. Further modification of the invention herein disclosed will occur to those skilled in the respective arts and all such modifications are deemed to be within the scope of the invention as defined by the appended claims. 

What is claimed is:
 1. A corneal examination method, comprising: (S1) reading an actual corneal deformation of a corneal of a tester from an ophthalmotonometer, wherein the ophthalmotonometer snaps corneal sectional images while blowing air onto the corneal, and generates the actual deformation according to the corneal sectional images; (S2) substituting a Young's modulus initial value and a damping coefficient initial value into a mathematical equation, wherein a calculated corneal deformation is obtained from a mathematical equation with the Young's modulus initial value and the damping coefficient initial value, and the mathematical equation is: ${w\left( {\theta,{\phi;t}} \right)} = {\sum\limits_{l = 1}^{\infty}\;{w_{l\; 0}\left( {\theta,{\phi;t}} \right)}}$ in which the term w_(l0) is the measured displacement on corneal surface ${w_{l\; 0}\left( {\theta,{\phi;t}} \right)} = {\frac{\left( {{2\; l} + 1} \right)}{2M_{l}\omega_{DI}}{{P_{l}\left( {\cos\;\theta} \right)} \cdot {\int_{0}^{\pi}{{g(\alpha)}{P_{l}(\alpha)}\sin\;\alpha\; d\;{\alpha \cdot {\int_{0}^{t}{e^{{- \xi_{l}}\omega_{Dl}t}\sin\;{\omega_{Dl}\left( {t - \tau} \right)}{s(t)}d\;\tau}}}}}}}$ where M_(l) is modal mass ${M_{l} = {{\rho \cdot h} + {R\left( {\frac{\rho^{+}}{l} + \frac{\rho^{-}}{l + 1}} \right)}}},$ modal frequency is ω_(Dl)=ω_(l)√{square root over (1−ξ_(l) ²)}, P_(l)(cos θ) is the Legendre polynomials, g(α) is the air pressure distribution on corneal surface, ξ_(l) is modal damping ratio, s(t) is the air puff time history, R is the corneal radius, h is the thickness of the diaphragm, ρ is the cornea density, ρ⁺ and ρ⁻ are densities inside and outside the cornea, l is the modal order, ${\omega_{l}^{2} = \frac{{\frac{\left( {E + T} \right) \cdot h^{3}}{12\left( {1 - v^{2}} \right)}\frac{{l^{3}\left( {l + 1} \right)}^{3}}{R^{4}}} + \frac{T \cdot h \cdot {l^{2}\left( {l + 1} \right)}^{2}}{R^{2}}}{{\left\lbrack {{\left( {l + 1} \right)\rho^{+}} + {l \cdot \rho^{-}}} \right\rbrack R} + {{l\left( {l + 1} \right)}{\rho \cdot h}}}},$ E is the Young's modulus, T is the tension on the diaphragm, and v is the Poisson's ratio, and images recorded from tonometer could provide the s(t), R, and h, and the densities, ρ, ρ⁺ and ρ⁻ are given from normal values, ξ_(l) is given 3˜5%, T is obtained from the intraocular pressure, v is assumed 0.4˜0.49; (S3) determining whether error amount between the calculated deformation and the actual deformation is a minimum or not; and (S4) determining the Young's modulus initial value and the damping coefficient initial value as a Young's modulus and a damping coefficient of the corneal, and judging whether the tester's corneal have disease, wherein the disease includes induced flaucoma, retinal lesion.
 2. The algorithm as claimed in claim 1, if the determining in the step (S3) is not, then back to the step (S2) and iterate in optimization until the error amount is a minimum.
 3. The algorithm as claimed in claim 1, wherein the minimum is 10⁻³.
 4. A corneal dynamic model system, comprising: an ophthalmotonometer, directing an external force to eyeball, and snapping corneal sectional images, generating an actual corneal deformation Dt; a calculating unit, connected to the ophthalmotonometer, wherein the calculating units acquires the actual corneal deformation from the ophthalmotonometer, and substituting a Young's modulus initial value and a damping coefficient initial value into a mathematical equation, wherein a calculated corneal deformation is obtained from a mathematical equation with the Young's modulus initial value and the damping coefficient initial value, and the mathematical equation is: ${w\left( {\theta,{\phi;t}} \right)} = {\sum\limits_{l = 1}^{\infty}\;{w_{l\; 0}\left( {\theta,{\phi;t}} \right)}}$ in which the term w_(l0) is the measured displacement on corneal surface ${w_{l\; 0}\left( {\theta,{\phi;t}} \right)} = {\frac{\left( {{2\; l} + 1} \right)}{2M_{l}\omega_{DI}}{{P_{l}\left( {\cos\;\theta} \right)} \cdot {\int_{0}^{\pi}{{g(\alpha)}{P_{l}(\alpha)}\sin\;\alpha\; d\;{\alpha \cdot {\int_{0}^{t}{e^{{- \xi_{l}}\omega_{Dl}t}\sin\;{\omega_{Dl}\left( {t - \tau} \right)}{s(t)}d\;\tau}}}}}}}$ where M_(l) is modal mass ${M_{l} = {{\rho \cdot h} + {R\left( {\frac{\rho^{+}}{l} + \frac{\rho^{-}}{l + 1}} \right)}}},$ modal frequency is ω_(Dl)=ω_(l)√{square root over (1−ξ_(l) ²)}, P_(l)(cos θ) is the Legendre polynomials, g(α) is the air pressure distribution on corneal surface, ξ_(l) is modal damping ratio, s(t) is the air puff time history, R is the corneal radius, h is the thickness of the diaphragm, ρ is the cornea density, ρ⁺ and ρ⁻ are densities inside and outside the cornea, l is the modal order, ${\omega_{l}^{2} = \frac{{\frac{\left( {E + T} \right) \cdot h^{3}}{12\left( {1 - v^{2}} \right)}\frac{{l^{3}\left( {l + 1} \right)}^{3}}{R^{4}}} + \frac{T \cdot h \cdot {l^{2}\left( {l + 1} \right)}^{2}}{R^{2}}}{{\left\lbrack {{\left( {l + 1} \right)\rho^{+}} + {l \cdot \rho^{-}}} \right\rbrack R} + {{l\left( {l + 1} \right)}{\rho \cdot h}}}},$ E is the Young's modulus, T is the tension on the diaphragm, and v is the Poisson's ratio, and images recorded from tonometer could provide the s(t), R, and h, and the densities, ρ, ρ⁺ and ρ⁻ are given from normal values, ξ_(l) is given 3˜5%, T is obtained from the intraocular pressure, v is assumed 0.45˜0.49, and the calculation unit determines whether error amount between the calculated deformation and the actual deformation is a minimum or not; and determines the Young's modulus initial value and the damping coefficient initial value as a Young's modulus and a damping coefficient of the corneal, and judging whether the tester's corneal have disease, wherein the disease includes induced glaucoma, retinal lesion.
 5. The system as claimed in claim 4, if the determining in the step (S3) is not, then back to the step (S2) and iterate in optimization until the error amount is a minimum.
 6. The system as claimed in claim 4, wherein the minimum is 10⁻³. 